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Abstract. Many mathematical models involve input parameters, which are not precisely 
known. Global sensitivity analysis aims to identify the parameters whose uncertainty has 
the largest impact on the variability of a quantity of interest (output of the model). One 
of the statistical tools used to quantify the influence of each input variable on the output is 
the Sobol sensitivity index. We consider the statistical estimation of this index from a finite 
sample of model outputs: we present two estimators and state a central limit theorem for 
each. We show that one of these estimators has an optimal asymptotic variance. We also 
generalize our results to the case where the true output is not observable, and is replaced 
by a noisy version. 

Resume. De nombreux modeles mathematiques font intervenir plusieurs parametres qui ne 
sont pas tons connus precisement. L'analyse de sensibilite globale se propose de selectionner 
les parametres d'entree dont I'incertitude a le plus d'impact sur la variabilite d'une quantite 
d'interet (sortie du modele). Un des outils statistiques pour quantifier I'influence de chacune 
des entrees sur la sortie est I'indice de sensibilite de Sobol. Nous considerons I'estimation 
statistique de cet indice a I'aide d'un nombre fini d'echantillons de sorties du modele: nous 
presentons deux estimateurs de cet indice et enongons un theoreme central limite pour chacun 
d'eux. Nous demontrons que I'un de ces deux estimateurs est optimal en terme de variance 
asymptotique. Nous generalisons egalement nos resultats au cas ou la vraie sortie du modele 
n'est pas observee, mais oil seule une version degradee (bruitee) de la sortie est disponible. 
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Introduction 



Many mathematical models encountered in applied sciences involve a large number of poorly-known 
parameters as inputs. It is important for the practitioner to assess the impact of this uncertainty on 
the model output. An aspect of this assessment is sensitivity analysis, which aims to identify the most 
sensitive parameters, that is, parameters having the largest influence of the output. In global stochastic 
sensitivity analysis (see for example (2^ and references therein) the input variables are assumed to 
be independent random variables. Their probability distributions account for the practitioner's belief 
about the input uncertainty. This turns the model output into a random variable, whose total variance 
can be split down into different partial variances (this is the so-called Hoeffding decomposition, see 
(3^ ) . Each of these partial variances measures the uncertainty on the output induced by each input 
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variable uncertainty. By considering the ratio of each partial variance to the total variance, we obtain 
a measure of importance for each input variable that is called the Sohol index or sensitivity index of 
the variable [27|; the most sensitive parameters can then be identified and ranked as the parameters 
with the largest Sobol indices. 

Once the Sobol indices have been defined, the question of their effective computation or estimation 
remains open. In practice, one has to estimate (in a statistical sense) those indices using a finite sample 
(of size typically in the order of hundreds of thousands) of evaluations of model outputs Indeed, 
many Monte Carlo or quasi Monte Carlo approaches have been developed by the experimental sciences 
and engineering communities. This includes the FAST methods (see for example Ql, [lH and references 
therein) and the Sobol pick-freeze (SPF) scheme (see (27l. [28|l. In SPF a Sobol index is viewed as the 
regression coefficient between the output of the model and its pick-freezed replication. This replication 
is obtained by holding the value of the variable of interest (frozen variable) and by sampling the other 
variables (picked variables). The sampled replications are then combined to produce an estimator 
of the Sobol index. In this paper we study very deeply this Monte Carlo method in the general 
framework where one or more variables can be frozen. This allows to define sensitivity indices with 
respect to a general random input living in a probability space (groups of variables, random vectors, 
random processes...). In this work, we study and compare two Sobol index estimators based on the 
SPF scheme; the first estimator, denoted by S"^, is well-known, the second, denoted by T§ has been 
introduced in For both estimators, we show convergence and give the rate of convergence; we 
also show that is optimal (in terms of asymptotic variance) amongst regular estimators which 
are functions of the pick-freezed replications - this feature is called asym pto tic efficiency and is a 
gen eralization of the notion of minimum variance unbiased estimator (see |32l | chapters 8 and 25 or 
\m for more details). 

The SPF method requires many (typically, around one thousand times the number of input variables) 
evaluations of the model output. In many interesting cases, an evaluation of the model output is 
made by a complex computer code (for instance, a numerical partial differential equation solving 
algorithm) whose running time is not negligible (typically in the order of a second or a minute) for 
one single evaluation. When thousands of such evaluations have to be made, one generally replaces 
the original exact model by a faster-to-run metamodel (also known in the literature as surrogate model 
or response surface [l|) which is an approximation of the true model. Well-known metamodels include 



Kriging [2j| , polynomial chaos expansion [SOj and reduced bases |12l . Il9l | , to name a few. When a 
metamodel is used, the estimated Sobol indices are tainted by a twofold error: sampling error, due to 
the replacement of the original, infinite population of all the possible inputs by a finite sample, and 
metamodel error, due to the replacement of the original model by an approximative metamodel. 

The goal of this paper is to study the asymptotic behavior of these two errors on Sobol index estimation 
in the double limit where the sample size goes to infinity and the metamodel converges to the true 
model. Some work has been done on the non-asymptotic error quantification in Sobol index estimation 
in earlier papers by means of confidence intervals which account for both samphng and 

metamodel errors. In this paper, we give necessary and sufficient conditions on the rate of convergence 
of the metamodel to the exact model for asymptotic normality of a natural Sobol index estimator to 
hold. The asymptotic normality allows us to produce asymptotic confidence intervals in order to assess 
the quality of our estimation. We also give sufficient conditions for a metamodel-based estimator to 
be asymptotically efficient. Asymptotic efficiency of an other Sobol index estimator has already been 
considered in [Hj. In this work, the authors were interested in the asymptotic efficiency for local 
polynomial estimates of Sobol indices. Our approach proposes an estimator which has a simpler form, 
is less computationally intensive and is more precise in practice. Moreover, we derive results also in 
the case where the full model is replaced by a metamodel. 

This paper is organized as follows: in the first section, we set up the notation, review the definition 
of Sobol indices and give two estimators of interest. In the second section, we prove asymptotic 
normality and asymptotic efficiency when the sample of outputs comes from the true model. These 
two properties are generalized in the third section where metamodel error is taken into account. The 
fourth section gives numerical illustrations on benchmark models and metamodels. 
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1. Definition and estimation of Sobol indices 

1.1. Exact model 

The output y e M is a functiou of independeut random input variables X G and Z e MP'^ . In 
other words, Y and {X, Z) are Hnked by the relation 

Y^f{X,Z) (1) 

where / is a deterministic function defined on T' C RPi+p^ , 'w^e denote by p = pi +p2 the total number 
of inputs of /. 

In the paper X' will denote an independent copy of X. We also write Y-^ = f{X, Z'). 

We assume that Y is square integrable and non deterministic (Var(y) ^ 0). We are interested in the 
following Sobol index: 

^ Var(E(y|X)) , , 

This index quantifies the influence of the X input on the output Y: a value of 5''''" that is close to 1 
indicates that X is highly influential on Y . 

Remark 1.1. All the results in this paper readily apply when X is multidimensional. In this case, 
S-^ is usually called the closed sensitivity index of X (see f2^]). 

Note that this separation between the input variables can be made without loss of generality, when one 
estimates Sobol indices independently. An ongoing work treats the case of joint Sobol index estimation. 

1.2. Estimation of 

The next lemma shows how to express S^ using covariances. This will lead to a natural estimator 
which has already been considered in Q. 

Lemma 1.2. Assume that the random variables Y and Y-^ are square integrable. Then 

Var(E(r|X)) = Cov(y,y^). 

In particular 

Remark 1.3. Using a classical regression result, we see that 

5^ = argmin |E((y-^-E(y^'^))- a(r-E(y)))H. (4) 

A first estimator. In view of Lemma [L^ we are now able to define a first natural estimator of 5*^ 
(all sums are taken for i from 1 to N): 

gX _ IV ^jYj^ ~ (iV S ^0 (iV S Yj^) 

where, for i = 1, . . . , A^: 

Y, = f[X,,Z,), Y.^' = f[X,,Z[), 

and {(^i, ^i)}i=i,- -.A^ ''^^'^ {{^ii Z[Y\i=i^...^N are two independent and identically distributed (i.i.d.) 
samples of the distribution of {X,Z), with {Zi}i independent of {Z'^}i. 

This estimator has been considered in where it has been shown to be a practically efficient esti- 
mator. 
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A second estimator. We can take into account the observation of {Fj"^}i<i<7v to make an 
estimation of E(y) and Var(y) which is expected to perform better than any other based on {yi}i<i<7v 
only. We propose the following estimator: 



N S 



Yi+Y.^ 



(6) 



This estimator has been introduced in We will clarify what we mean when saying that 

performs better than S§ in Proposition 12. 3( Section [2?2l and Subsection HJ] 

Remark 1.4. Note that the empirical variances in and can be rewritten as: 



S 



N 



(7) 
(8) 



where: 



Y + Y^ 



TV ^ ' N 2 

The use of these formulae enables greater numerical stability (ie., less error due to round-off s). The 
Kahan compensated summation algorithm fl2] may also be used on these sums. However, we will 
use definitions ([5]) and ^ for the mathematical analysis of and T^ . This analysis is of course 
independent of the way the estimators are numerically computed in practice. 



2. Asymptotic properties: exact model 

2.1. Consistency and asymptotic normality 

Throughout all the paper, we denote by A/fc(jU, S) the fc-dimensional Gaussian distribution with mean 
yU and covariance matrix E, and, given any sequence of random variables {i?„}„gN, we note 

N 



ri=l 



Proposition 2.1 (Consistency). We have: 



qX a-s-^ qX 

D AT ? D 



7V- 



T, 



X g-^^ gx 

N~ 



N 



(9) 

(10) 



Proof. The result is a straightforward application of the strong law of large numbers and that E,{Y) = 
E(y^) and Var(y) = Var(y^). □ 

Proposition 2.2 (Asymptotic normality). Assume that E{Y*) < 00. Then 



and 



whe 



(11) 
(12) 



Var {{Y - E{Y)) [{Y^ - E{Y)) - (Y - E(y))]) 
(Var(r))^ 
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2 _ Var {{Y - E{Y)){Y^ - E{Y)) - /2 {{Y - E{Y))^ + {Y^ - E(r))')) 
~ (Var(r))2 • 

Proposition 2.3. The asymptotic variance of is always less than or equal to the asymptotic 
variance of , with equality if and only if 5^ ~ or 5^ ~ 1 . 



To prove this Proposition, we need tire following immediate Lemma: 

Lemma 2.4. Y andY^ are exchangeable random variables, ie. (Y,Y-^) = {Y^,Y). 

2.2. Asymptotic efficiency 

In this section we study the asymptotic efficiency of S§ and . This notion (see [s^l, Section 25 for 
its definition) extends the notioir of Cramer-Rao bound to the semiparametric setting and enables to 
define a criteria of optimality for estimators, called asymptotic efficiency. 

Let V be the set of all cumulative distribution functions (cdf) of exchangeable random vectors in 
i^(]R^). It is clear that the cdf Q of a random vector of L^(R^) is in V if and only if Q is symmetric: 

Q{a,b) ^ Q{b,a) y{a,b)eR'^. 
Let P be the cdf of {Y, Y^). We have P eV thanks to Lemma [221 

Proposition 2.5 (Asymptotic efficiency). {T^}n is asymptotically efficient for estimating 5^ for 

PeV. 

We will use the following Lemma, which is also of interest in its own right: 

Lemma 2.6 (Asymptotic efficiency in V). (1) Let $i : M ^ R &e a function in L^{P). The 
sequence of estimators given by: 



1 ^ <i>i(yo + $i(y,^) 



is asymptotically efficient for estimating E{^i(Y)) for P Cz V . 
(2) Let $2 : — > M 6e a symmetric function in L'^{P). The sequence {^%} given by: 



N 

is asymptotically efficient for estimating E{^2(Y,Y'^)) for P (zV . 

3. Asymptotic properties: metamodel 

3.1. Metamodel-based estimation 

As said in the introduction, we often are in a situation where the exact output / is too costly to be 
evaluated numerically (thus, Y and Y^ are not observable variables in our estimation problem) and 
has to be replaced by a metamodel /, which is a faster to evaluate approximation of /. We view this 
approximation as a perturbation of the exact model by some function 5: 

Y = f{X,Z)^f{X,Z) + 6, 

where the perturbation S = 5(X,Z,(^) is also a function of a random variable f independent from X 
and Z . 

We also define, as before 

f^ = /(X,Z'). 

Assuming again that Y is non deterministic and in L^, we can consider the following Sobol index, 
with respect to the metamodel: 

~, ^ Var(E(yjX)) 
Var(r) 
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and its estimators: 



N 



(14) 



N S 



Y^^ + {Y.^ 



Yi+Y.^ 



(15) 



The goal of this section is to give sufficient conditions on the perturbation S for and to satisfy 
asymptotic normality (Subsection I3.2p . and to be asymptotically efficient (Subsection I3.3p . with 
respect to the Sobol index of the true model . 



3.2. Consistency and asymptotic normality 



In the first Subsection (|3.2.ip we suppose that the error term S does not depend on A''. In this case, if 
the Sobol index of the exact model is different from the Sobol index of the metamodel, then neither 
consistency nor asymptotic normality are possible. In the second subsection (|3.2.2p . we let S depend 
on N and we give conditions for consistency and asymptotic normality to hold. 



3.2.1. First case : S does not depend on N 



Remark 3.1. If 

Indeed, we have 



S 7^ then neither Sj^ nor are consistent for estimating S 



X 



-,x 



X 



-yX\ 



The first term converges to almost surely by Proposition 
nonzero by assumption. The same holds for T^ . 



applied to S^. However, the second is 



This remark shows that a naive consideration of the metamodel error (ie., with fixed metamodel) 
is not satisfactory for an asymptotic justification of the use of a metamodel. More specifically, it is 
impossible to have asymptotic normality for Sj^ and in any nontrivial case if S does not vanish 
(in some sense) asymptotically. This justifies the consideration of cases where 6 depends on A^, and 
this is the object of the next subsection. 



3.2.2. Second case : Var Sn converges to as N 



oo 



We now assume that the perturbation i5 is a function of the sample size N. This entails that /, as 
well as Y, Y^ and depend on N. We emphasize this dependence by using the notations Sn, fx, 
Yx, ^N ■ We keep, however, using the notations S-§ and for the estimators of S-^ defined at (fM)) 
and (fTSt. 



Assumption. We suppose that /jv 
Proposition 3.2. We have 



f 



Sn — > c for some constant c. 

Ar-s.+oo 



W-S.+00 

Proposition 3.3. Assume there exist s > and C > such that 



ViV, E 



Yn 



4+s 



< c. 



(16) 



The 



N S 



^x 



^x 

X 



X 



c ^ 

TV-s-oo 

x\ 



Ml {0,4) 
Ml (0,4) 



(17) 
(18) 



where ag and are the asymptotic variances of and given, respectively, in (jll[) and (fT2 
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We are actually interested in the asymptotic distribution of \fN yS-§ — S-^j . In the remaining of the 
subsection, we will show that this convergence depends on the rate of convergence to of Var((5Ar). 

Theorem 3.4. Let: 

CsM = 2Var(y)i/2 [CoTi{Y,d§) - Corr(y, y^)Corr(r, Jat)] +Var(5Ar)i/2 [CoTT{dN,S^) - Corr(y,y^)] , 
for 6^ — 6n{X,Z'), and, given any random variables A and B of nonzero variance: 

Cov(^, B) 



Corr(A, B) 



WarAVarS' 



Assume that Cs.n does not converge to 0. 

(1) If Var{S]y) =o(-^), then asymptotic normalities of and for S'^ hold, i.e. 



mSj^-S'') Af{0,al) (19) 

and: 

VN{f^-S'') AA(0,4). (20) 

A^— >-+oo 

(2) IfNY&iidN) oo, then (HH) and 

(3) If Cs,N converges to a nonzero constant C and 7 G M so that Var((5jv) = + (■^), then: 



and: 

ViV(f#-^^) —> AA(7,4)- 

A*— >-+oo 

Remark 3.5. Obviously, if Cs.n converges to 0, then asymptotic normalities of S*^ and hold 
under weaker assumptions on Var((5jv)- 

3.3. Asymptotic efficiency 

Proposition 3.6 (Asymptotic efficiency for the metamodel). Assume 
(1) 3s > 0, C > s.t. yN, E (\Y\'^^') <C andE (\Y\^^') < C ; 



(2) NYariSN) ~^ ; 

(3) VnE{6n)^0. 

Then |'^jv | asymptotically efficient for estimating S'^ . 

Remark 3.7. By Minkowski inequality, the first hypothesis implies E((5^'*) < 2C^ and the asymp- 
totic normality by Lemma \cl.3\ and Theorem 



4. Numerical illustrations 

In this section, we illustrate the asymptotic results of Sections 12.11 and 13.21 when the exact model is 
the Ishigami function 

f{Xi,X2,X3) =sinAi + 7sin2A2 + 0.1A|sinAi (21) 

for (Aj)j=i^2,3 are i.i.d. uniform random variables in [— 7r;7r]. In this case, all the integrability 
conditions are satisfied (we even have Y <E L°°). 

The Sobol index of / with respect to input variable Xi is defined in 1^ ior X = Xi and Z = 
{X2 , A3) ; we denote it by . Similarly, S'^ (resp. S^) is 5"^ obtained taking A = X2 and Z = (Ai , A3) 
(resp. A = A3 and Z = (Ai, A2)). 
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Figure 1. Empirical coverages of asymptotic confidence intervals for (left), 5^ 
(center) and (right), as a function of the sample size. The estimator is used. 




Figure 2. Empirical coverages of asymptotic confidence intervals for (left), 
(center) and (right), as a function of the sample size (for the exact model). The 
Tn estimator is used. 



Exact values of these indices are analytically known: 

= 0.3139, = 0.4424, = 0. 

For a sample size N, a risk level a e]0; 1[ and for each input variable, a confidence interval for S'^ 
{S^ being one of S^, S'^ or S^) of confidence level 1 — a can be estimated - using evaluations of the 
true model / - by approximating the distribution of (or ) by its Gaussian distribution given 
in Proposition 1111 using empirical estimators of the asymptotic variances stated in this Proposition. 

In the case where only a perturbated model (metamodel) /at = / + (5jv is available, a confidence 
interval can still be estimated by using the (or T^) estimator. 

Thanks to Proposition 13.31 the level of the resulting confidence interval should be close to 1 — q for 
sufficiently large values of N if (and only if) YarSN decreases sufficiently quickly with N. 
The levels of the obtained confidence interval can be estimated by computing a large number R 
of confidence interval replicates, and by considering the empirical coverage, that is, the proportion 
of intervals containing the true index value; it is well known that this empirical coverage strongly 
converges to the level of the interval as R goes to infinity. 

In the next subsections, we present the estimations of the levels of the confidence interval for the 
Ishigami model (PT|) using the true model (Subsection 14. ip . and, with various synthetic model per- 
turbations (Subsections 14.21 and 14. 3p . as well as RKHS (Kriging) metamodels (Subsection 14. 4p and 
nonparametric regression metamodels (Subsection 14. 5p . We begin by comparing Sn and Tjv on the 
exact model (Subsection 14. ip . then we illustrate the generalization to the metamodel case on the 
widespread estimator Sn\ the condition to ensure asymptotic normality in the metamodel is the same 
for Sn and Tjv. All simulations have been made with R = 1000 and a ~ 0.05. 



4.1. Exact model 

Figure [T] shows the empirical coverage of the asymptotic confidence interval built using the 5^ esti- 
mator, plotted as a function of the sample size N. The theoretical level 0.95 is represented with a 
dotted line. Figure [2] does the same using the estimator. 
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Figure 3. Lengths (rescaled by vN) of the estmiated 95% confidence intervals for 
(left), S'^ (center) and (right), as functions of the sample size (for the exact 
model). In solid line: length of the interval built from T/v estimator; in dotted line: 
length of the interval built from Sn estimator. 
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Figure 4. Empirical coverages of the asymptotic confidence intervals for S*^, 5^ and 
S^, as a function of /3 (for the Gaussian-perturbated model). 



We see that the coverages get closer to the target level 0.95 as N increases, thereby assessing the 
reliability of the asymptotic confidence interval. 

Figure [3] compares the efficiency of and by plotting the confidence interval lengths for the two 
estimators, as functions of the sample size. As the lengths for both estimators are 0{1/VN), we plot 
the lengths multiplied by \/N. We see that always produce smaller confidence intervals, except 
for X3 where the lengths are sensibly the same; this conclusion fully agrees with Proposition 12.31 



4.2. Gaussian-perturbated model 

We consider a perturbation of the original output /: 

5^ 



where /? > and ^ is a standard Gaussian. 

The perturbation (5jv = 5 ^^^^2 leads to Var^Ar oc N"^ . Since: 



C, = 0(var(M'/')=0(^-^/'), 



the proof of Theorem l3.4l shows that 5jv is asymptotically normal for 5 if /? > 1/2. For indices relative 
to Xi and X2, this sufficient condition is also necessary, as Cs is actually equivalent to iV~'^/^. For 
X3, we have = so that Sn is asymptotically normal for S for any positive j3. 

This is illustrated for N ~ 50000 in Figure [H We see that the empirical coverages of the confidence 
interval for and jump to 0.95 near /3 = 1/2, while, for , this coverage is always close to 0.95. 
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Figure 5. Empirical coverages of the asymptotic confidence intervals for S^, and 
S^, as a function of f3 (for the Weibull-perturbated model). 



4.3. Weibull-perturbated model 

We now take a different perturbation of the output: 



hWXl 

Ar-9/2 



where W is Weibull-distributed with scale parameter A = 1 and shape parameter k = 1/2. Here, 
the perturbation depends on the inputs and, as for every input variable, Cg.N does not converge to 
zero, Theorem 13.41 states in particular that Sn is asymptotically normal for S for /3 > 1. Again, this 
property is suggested for N = 50000 by the plot in Figure [5] 

4.4. RKHS metamodel 

In this part, we discuss the use of a reproducing kernel Hilbert space (RKHS) interpolator (23. [25l. [26t 
as metamodel /. Such metamodels (also known as Kriging, or Gaussian process metamodels) are 
widely used when performing sensitivity analysis of time-expensive computer codes (l6| . Note that, 
according to Q, analytical formulae are in some cases (e.g., uniform or gaussian distributions for 
the inputs) available for Sobol indices computation, avoiding the necessity to use a Monte-Carlo 
scheme. In this paper, we chose to perform Monte-Carlo estimation on an RKHS metamodel so 
as to illustrate our theoretical results. Moreover, the Monte-Carlo approach is more flexible and 
can be applied for complex inputs' distributions. The interpolator depends on a learning sample 
{(di, /(di)), . . . , (dn, /((i„))}, where the design points T) = C V are generally chosen 

according to a space-filling design, for instance the so-called maximin LHS (latin hypercube sampling) 
designs. Increasing the learning sample size n will increase the necessary number of evaluations of 
the true model / (each evaluation being potentially very computationally demanding) to build the 
learning sample, but will also enhance the quality of the interpolation (i.e. reduce metamodel error). 

The error analysis of the RKHS method [l5, 25| shows that there exist positive constants C and /C, 
depending on /. so that: 



f{u) - f{u) 



where: 



sup mm I 



for a given norm 1 1 • 1 1 on T-" . 

The quantity hxi;p can be linked to the number of points n*{e) in an optimal covering of T): 

n*{e) = min{p e N* \ 3{di, . . . ,dp) e V s.t. Vu e P, 3i e {1, . . . ,p} satisfying \\u - di\\ < e}. 

In other words, n*{e), known as the covering number ofV, is the smallest size of a design V satisfying 
hv,'P < £• 

It is known that, when T' is a compact subset of (in our context, p = pi + P2 is the number of 
input parameters), there exist constants A and B so that: 



Ae-P < n*{e) < Be'P. 
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Figure 6. Estimation of the Kriging metamodel error variance (log. scale) as func- 
tion of the learning sample size n. 



Hence, assuming that an optimal design of size n is chosen, we have, for a constant B' : 
and we have the following pointwise metamodel error bound, for constants C and K' : 



fin) - f{u) 



<Ce 



which obviously leads to an integrated error bound on the variance of the metamodel error: 



VarJ < Ce 



for suitable constants C and k. 



Numerical illustration 



We illustrate the properties of the RKHS-based sensitivity analysis using the Ishigami function ([21]) 
as true model, maximin LHSes for design points selection. RKHS interpolation also depends on the 
choice of a kernel, which we choose Gaussian all the way through. All simulations have been made 
with the R software [23, together with the Ihs package Q for design sampling and the mlegp package 
for Kriging. 

Figure [HI which shows an estimation (based on a sample of 1000 metamodel errors) of the (logarithm 
of) variance of metamodel error, plotted against the cubed root of the learning sample size n^^^. Using 
an exponential regression, we find that: 



Var((5) « Ce" 



(22) 



where: 

k = 1.91 

Now, if we let the learning sample size n depend on the Monte-Carlo sample size N by the relation: 

n = (alnA^)^ 
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a 


TV 


n 


Coverage for 5*^ 


Gov. for S'^ 


Gov. for S'-^ 


A 


3000 


33 


0.1 





0.7 


A 


4000 


37 


0.08 





0.78 


A 


6000 


43 


0.26 


0.3 


0.88 


A 


10000 


51 


0.28 


0.18 


0.78 


A 


20000 


77 


0.28 


0.1 


0.59 


.6 


3000 


111 


0.79 


0.37 


0.9 


.6 


4000 


124 


0.8 


0.7 


0.94 


.6 


10000 


169 


0.92 


0.82 


0.94 


.6 


20000 


210 


0.93 


0.85 


0.95 


.7 


3000 


177 


0.93 


0.88 


0.93 


.7 


4000 


196 


0.9 


0.91 


0.94 


.7 


6000 


226 


0.94 


0.93 


0.97 


.8 


4000 


293 


0.95 


0.95 


0.95 



Table 1 . Estimation of the asymptotic coverages for the RKHS Ishigami metamodel. 
Empirical coverages are obtained using 100 confidence interval replicates. Theoretical 
coverage is 0.95. 



for a > 0, Theorem 13.41 suggests that the metamodel-based estimators of the sensitivity indices are 

when N — > +00, that is a > i, or 

a > 0.52, (23) 



asymptotically normal if and only if N — > when N — > +00, that is a > i, or 



according to our numerical value for k. 

Even if it has not been rigorously proved that this condition is necessary and sufficient (due to the 
estimation of k and the fact that (|22p provably holds, possibly with different constants, as an upper 
bound), one should observe in practice that the behavior of the empirical confidence intervals for large 
values of N changes as this critical value of a is crossed. Table [1] below shows the results obtained for 
different subcritical and supercritical values of a (i.e., (|23p does not hold, or hold, respectively), and 
provides a clear illustration of this fact. 

4.5. Nonparametric regression 

In this section, we consider the case where the true model / is not directly observable, but is only 
available through a finite set of noisy realisations of: 

/noisy(A) = /(A) + Ci, i = l,...,n 

where T) = {Di = {Xi, Zi))^^^ ^ are independent copies of {X, Z), and {ei}i=i....,n are independent, 
identically distributed centered random variables. 

As discussed in Section [3.2.11 one should expect that the Sobol index estimator computed on /noisy 
are not asymptotically normal for the estimation of the Sobol indices of / (as Var(ei) is fixed). This 
motivates the use of a smoothed estimate of /, which we will take as our perturbated model / = /p. 
We consider the Nadaraya- Watson estimator: 

^-^ff:(V^')^7^^-^ fx:i^..(^^-A)^o 

else. 

where Kh is a smoothing kernel of window h G R^; for instance Kh is a Gaussian kernel: 

K,iv) = e.J-j2^^] (24) 



where the norm ||-|| is the Euclidean norm on MP. 
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It is known that, under regularity conditions on /, and a 7i-dependent appropriate choice of /i, the 
mean integrated square error (MISE) of / satisfies: 

^(^[f{u)-fv{u)y^du<C'n-\ (25) 

for a positive constant C" and a positive 7 (which depends only on the dimension p and the regularity 
of /), and where E-p denotes expectation with respect to the random "design" V. 

Now, by Fubini-Tonelli's theorem, we have: 

Ev(^[f{u)-fv{u)y^du = Ev(^l [f{u)-fv{u)yduy (26) 



Bv using ([25 ]) . (P5|) and applying Markov's inequality to the positive random variable J ^/(u) — /x)(u)^ du, 
we have that, for any e > 0, 

T^/ y (/(")-/p(ti))'cl«< > 1-e. 

Hence, for a fixed risk e > 0, there exist C > and 7 > so that: 

'fviu)-f{u)Ydu<Cn-^ (27) 



holds with probability greater than 1 — e (with respect to the choice of V) . 

We recall that the quantity we have to consider in order to study asymptotic normality of Sobol index 
estimator on the metamodel is: 

Var(5) = J (^f{u)-fv{u)y du- (^J (^f{u)-fv{u)) du^ 

and that, obviously, 

Var(5) < I (^f{u)-Mu)ydu. 

This gives, by making use of (P7)) : 

Var (S) < Cn-'^ (28) 

with probability greater than 1 — e. 

In most cases of application, the design V is fixed. In view of (|28p . it is reasonable to suppose that 
there exist C > and /3 > so that: 

Var [S) < Cn'^ 
and we make n depend on N by the following relation: 

n = N", 

for a > 0. By Theorem 13.41 the estimator sequence {S'Ar} is asymptotically normal provided that 
NYa.i-{6N) 0, that is: a > -g. 

Numerical illustration 

We now illustrate this property using the Ishigami function (|2ip as true model, and a Gaussian white 
noise of standard deviation 0.3 (yielding to a signal-to- noise ratio of 90%). 

The nonparametric regressions are carried using a Gaussian kernel p4| , the R package np Q , together 



with the extrapolation method of [2l| for window selection and the FIGtree [18| C++ library for 
efficient Nadaraya- Watson evaluation based on fast gaussian transform. 

Figure [3 which shows an estimation (based on a test sample of size 3000) of Var{5) in function of n, 
and a power regression shows that: 

Var {S) w Cn"^ 
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Figure 7. Estimation of the nonparametric regression error variance (log. scale) as 
function of the learning sample size n (Subsection 14 . 5 p . 



a 


N 


n 


Coverage for 


Gov. for 


Gov. for S'-" 


0.8 


1000 


252 


0.25 


0.01 


0.94 


0.8 


2000 


438 


0.05 


0.02 


0.86 


1.1 


1000 


1996 


0.95 


0.97 


0.96 


1.1 


2000 


4277 


0.95 


0.93 


0.96 


1.2 


1000 


3982 


0.93 


0.95 


0.96 


1.2 


2000 


9147 


0.96 


0.97 


0.95 


1.3 


1000 


7944 


0.95 


0.99 


0.94 


1.3 


2000 


19559 


0.95 


0.95 


0.96 



Table 2. Estimation of the asymptotic coverages for the Ishigami nonparametric 
regression. Empirical coverages are obtained using 100 confidence interval replicates. 
Theoretical coverage is 0.95. 



with (3 = 0.86. This gives an estimate of 1.16 as the critical a for asymptotic normality. 

As in the RKHS case, we performed estimations of the coverages of the asymptotic confidence interval 
for several values of a and N; the results are gathered in Table [5J We see that, first, the condition 
a > 1.16 implies correct coverages, and, second, the condition also seems to be near- necessary to have 
asymptotic normality. We also remark that, for the asymptotic normality to hold, the necessary num- 
ber of noisy model evaluations is asymptotically comparable to the Monte-Garlo sample size (while, 
in the RKHS case, the necessary number of true model evaluations was asymptotically negligible with 
respect to the Monte-Garlo sample size): this shows that the nonparametric regression is suitable 
in the case of noisy but abundant model evaluations, while RKHS interpolation is clearly preferable 
when the true model output is costly to evaluate (i.e. few model outputs are available). 



5. Appendix: Proofs 

Proof of Lemma \1.2\ On one hand, since Y = Y^ (that is, Y and Y^ have the same distribution), 
we have 

Gov(y, Y^) = E(yy-^) - ^{Y)^{Y^) = "EiYY^) - ^{Yf. 
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On the other hand, Y and Y-^ are independent conditionally on X, so that 

E{YY^) = E{E{YY^\X)) = E{E{Y\X)E{Y^\X)) = E{E{Y\Xf). 

□ 

Proposition \2.S[ ' Proof of (jlip . We begin by noticing that is invariant by any centering (transla- 
tion) of the Yi and Y/^ . To simplify the next calculations, we suppose that they have been recentred 
by -E(r). By setting: 

= {{Y,-E{Y)){Y/' -E{Y)), Y,-E{Y), Y/" -E{Y), {Y,-E{Y))Y, (29) 
this implies that: 

S§ = ^s{Un) 

with: 

, . X — yz 

'$s[x,y,z,t) = J 

t - 

The central limit theorem gives that: 

N^OQ 



where F is the covariance matrix of Ui and: 



/Cov(y,y^)\ 




V Var(y) J 



NiS^-S"") A/'i(0,.g^rg) 

N—^oo 



The so-called Delta method [32| (Theorem 3.1) gives: 
where: 

<7 = V*s(m)- 

Note that since by assumption Var(F) 7^ 0, ^5 is differentiable at fi and we will see that g^Tg ^ 0, so 
that the application of the Delta method is justified. By differentiation, we get that, for any a;, y, z, t 
so that t ^ y"^: 



V^s{x,y,z,t) 
so that, by using ([3]): 



1 -z{t - y"^) + {x-yz)-2y y x ~ yz ^ ^ 



t-y^' {t~y^Y ' t-y2' - y^f 



0, 0, * 



Var(y) ' ' ' Var(y) J 



Hence 



Var((r-E(r))(y^ -E(r))) (s^Y 



(Var(y))2 
1 



Gov ((y - E{Y)){Y^ ~ E(r)), (y - E{Y) f) 



Var ((y - E(y))(y-^ - E(y))) + Var {S^ ((y - E{Y) f)) 



(Var(y))2 ^ 

-2Cov ((y - E(y))(y'^ - E(y)), s^{y - E{Y)f) 
Var ((y - E(y)) [(y^ - E(y)) - ^^(y - E(y))] ) 

(Var(y))2 ' 
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which is the announced result. 

Proof of (fT2|) . As in the previous point, it is easy to check that is invariant with respect to 
translations of Y, and Y/^ by -E{Y). Thus, = 1- (Wn) with: 



X - {y/2f 



z/2-(y/2)2 
and: 

= ((y, - E(F))(K,^ - E(F)), (y, - E(r)) + (k,^ - E(r)), (k, - E{Y)f + (y,^ - ny)ff ■ 

(30) 

The result follows from the delta method. 

□ 

Proof of Proposition \2.Si We have that the expressions in (jTI]) and of cr| and (tI^ are translation- 
invariant, so that we assume without loss of generality that E(y) = 0. By expanding the variances 
and using the exchangeability of Y and Y^ ^ we have that (Var(y))^ (cr| — is equal to: 

(Var(y))2 {al - 4) ^ (Var(y2) - Gov (r^ (r^')^)) . 



We now use Cauchy-Schwarz inequality to see that: 

Gov {Y"^, {Y^ f) < y^Var(r2)Var(Y^)2) = Var (F^) 

so the second term is always non-negative. This proves that the asymptotic variance of is greater 
than the asymptotic variance of . 

For the equality case, we notice that S'^ = implies the equality of the asymptotic variances. If 
S-^ 0, equality holds if and only if there is equality in Cauchy-Schwarz, ie. there exists fc € R so 
that: 

= k{Y^f almost surely 

by taking expectations and using Var(F) = Var(y'^) we see that fc = 1 necessarily, hence Y = Y^ 
almost surely, and S-^ = 1 thanks to ([3|). □ 

Proof of Lemma lKR Let, for g G L?{P) and t G R, P^ be the cdf satisfying: 

dP/ = (1 + tg)AP. 
It is clear that the tangent set of P at P is the closure of: 

Vp = {g bounded, ¥.{g{Y,Y^)) = and g{a,h) = g{b,a) V(a, fe) e M^}. 



Let, for Q eV: 

*i(g) =EQ($i(y)) and ■^2{Q)=Eq{^2{Y,Y'')). 

We recall that Eg denotes the expectation obtained by assuming that the random vector (Y^Y-^) 
follows the Q distribution. 

Following [s^l Section 25.3, we compute the efficient influence functions of and ^2 with respect to 
V and the tangent set Vp. These empirical influence functions are related to the minimal asymptotic 
variance of a regular estimator sequence whose observations lie in V (op.cit.. Theorems 25.20 and 
25.21). Let g e Vp. 

(1) We have 



«'i(Pf)-*i(P) 



= E 



Ep {<^i{Y)g{Y,Y'')) 

$i(r) + $i(r^) 



E(<I>i(y)))g(y,F^) 
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As: 



c:;^5i2>±£ifi:i)_E(^..(y))cPp. 



it is the efficient influence function of ^'i at P. Hence tlie efficient asymptotic variance is: 



Var($i(r) + 



As, by tlie central limit theorem, {^'jv} clearly achieves this efficient asymptotic variance, it 
is an asymptotically efficient estimator of \E'i(P). 
(2) We have: 

^pn-^p) ^ Ep(a>2(y,y^).g(y,F^)) 

Thanks to the symmetry of $2, we have that 

belongs to Vp, hence it is the efficient influence function of 4*2 • So the efficient asymptotic 
variance is: 

Ep(^($;S)') =V&r{^2{Y,Y'')) , 
and this variance is achieved by { $^ } . □ 
Proof of Provosition \2.5l By Lemma 12.61 we get that: 



\ N ^ ' ' N ^ 2 N ^ 

\ i=l 1=1 1=1 



is asymptotically efficient, componentwise, for estimating 

U = (E{YY^), E(y), E(r2)) (32) 

in V. 

Using Theorem 25.50 (efficiency in product space) of (32| . we can deduce joint efficiency from this 
componentwise efficiency. 

Now, let be the function defined by: 

U'(x,y,z) = 2 

and 'f is difFerentiable on: 

R''\{{x,y,z)\z^y^}, 

Theorem 25.47 (efficiency and Delta method) of [s^l implies that {"^ (Un)} is asymptotically efficient 
for estimating "^{U) for P eV. The conclusion follows, as '^{Un) = and *([/) = . □ 

Proof of Proposition \3.Si We clearly have that Yn — > Y + c. 

N-^+oo 

We deduce that: 



and 



Var (Yn] — ^ Var(y + c) = Var(y) 

V / Af->+oo 



E{Yn\Z) — > E{Y\Z) + c in 

>+oo 
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From this last convergence we get 

YnT (e{Yn\Z)] — > Var(E(y|Z)). 

This proves that = Var (e{Yn\Z)^ /Var (Yn^ converges to = Var (E(y|Z)) /Var(r) when N 
goes to +00. □ 

Proof of Proposition\KE Proof of Let 

UN,^ - (^{YN,^ - E(y))(r#, - E{Y)), Fat., - E{Y), f^, - E{Y), (Yn., - E(y)) 

and 

AT 



i=l 

Using the Lindeberg-Feller central limit theorem (see e.g. [s^ 2.27, with Y^^i = UN.i/VN), we get: 



where F is the covariance matrix of the Ui vector defined in ([29]) . 

The use of this central limit theorem is justified by the fact that, under assumption (jl6p of uniform 
boundedness of moments of Y/v, there are s' > and C such that: 

V7V, E{\\UnA?^"')<C' 

where ||-|| is the standard Euclidean norm. 
This ensures 

w, ^ n 1T7/Il7~r Il2, _ 



Ve>0, E{\\UNA?\fj.A\>.^)^^- 



Then 



This shows that for each 



\UnAY 

is uniformly integrable, hence, the variance-covariance matrix 



N 

of [/^r i converges to F when N — > +00. As Un i — > Ui, the same convergence holds in L and the 

covariance matrices of Una converge (as N — > +00) to F, the covariance matrix of Ui. 

We conclude the proof by applying the Delta method as for the exact model (cf the proof of Propo- 
sition d^l). 

Proof of llH]). We set: 



WN,^ = [{Y, - E{Y))(Y..^ - E(y)), (F, - E{Y)) + (F,^ - E(y)), (F, - E{Y)f + (F,^ - E{Y)f) ' 
As in the previous point, the Lindeberg-Feller theorem can be applied to |vF/v^i| to yield the conver- 



gence: 



N {Wn -"^(Wi.i)) ^ AA3(0,I]) 



where S is the covariance matrix of Wi defined in ([5U|) . The conclusion follows again by an application 
of the Delta method as in the proof of Proposition 12.21 □ 



Proof of Theorem \3.4\ The following decompositions: 

Vn{S§-S^) = Vn{S§~S^) + Vn{S^ -S^) (33) 
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X qX\ _ qX\ I ^fTjiaX qX 



N—>--\-oo 



'N{T^ -S^) = VN{T^ -S^) + VN{S^ ~ ) (34) 
make obvious that if \/N{S^ - S^) goes to some constant k then 



and: 

The second point of the theorem is now clear from the proof of Proposition 13.11 

The remaining of the theorem is an immediate consequence of Lemma 15.11 below. □ 
Lemma 5.1. We have: 

r- O ((iVVar(5^))^/^ 



Var(y) + o(l) 



Proof of l5.il We have: 



cX 



Cov(rjv,F#) Cov(r,y^) 



VarYiv Var(y) 



Cov(F, Y^) + 2Cov(r, 5^) + Cov((5a,, 5^) Cov(r, Y^) 



Var(y) + 2Cov(r,(5Ar) + Var((5jv) Var(y) 
Var(y) (2Cov(r, S§) + Cov{Sn,5^)) - Cov{Y, Y^) (2Cov(r, Sn) + Var(<5jv)) 
Var(r) (Var(y) + 2Cov(r, Sn) + Var(,57v)) 
Var(<SAr)i/2C5„ 



Var(y) + 2Cov(r, Sn) + Var((5Ar) 

and: 

Var(y) + 2Cov(y, 5n) + Var((5Ar) = Var(y) + o(l). 
Finally, Cs^n is uniformly bounded because Va.r{SN) goes to and Var(y) is a constant. □ 

Proof of Proposition \3.6l We will use the following lemma. 

Lemma 5.2. For all N G N*, let {ZN,i)i=i,...,N be a sequence of i.i.d variables such that 



(1) VNE{ZNa) 0; 

(2) Var(Zjv,) — > 0. 

Then 

-^YZn. a 0. 

^ 2—1 

The lemma follows after the following decomposition: 

1 ^ / 1 ^ \ 

^ i—l \ i—1 / 

Let Un and U be the vectors defined in the proof of Proposition 12.51 in (PT|) and ([5^ , respectively, 
and: 

\ i=l i=l i=l / 



We will show that: 



(Un - Un) 0. (35) 
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By Theorem 25.23 of [13 and the fact that (f/jv) is asymptotically efficient for U (shown in the proof 
of Proposition 12. 5p . this implies that (Un^ is asymptotically efficient for U, and the end of the proof 
of Proposition 12.51 shows the annomiced result. 

To prove ([35|) . it is sufficient to prove componentwise convergence. We will treat the second and the 
third components, as the result holds in the same way for the other. 

For the second component, we have 

1 ^ ~ 1 " 

goes to (in probability) by the previous lemma. The same holds for '^iLiO^N i ~ ^i^)- 
For the third component, we have 

N N N 

Now by assumption, 

and by Cauchy-Schwarz inequality, 

V'AY{SN.Yi)^E{6%^,Y^) ~ {E{SNaY.)f < ^E{6%.^)E{Y^) + E(4.,)E(y/) < CE{S%)^^^. 

p /I . . . 

By assumption, for all i, 5N,i — > 0. Hence, the same convergence holds about ^. Since 5^ is in 

L^+^. then {S%} is uniformly integrable and we get the convergence of E{5^) to when +oo. 

We conclude by the lemma above. Again, the same convergence occurs for X^i^i ((^iv"i)^ ~ (Yi^)'^)- 

□ 

Conclusion 

We have shown that the Sobol index estimator considered in this paper is asymptotically normal and 
asymptotically efficient. We also proved that these two properties are robust with respect to the use 
of a perturbated model, provided that the perturbation has a variance which decays fast enough. The 
asymptotic normality property can be used to produce approximate confidence intervals for the Sobol 
indices; we have presented numerical experiments asserting the reliability of these confidence intervals. 

This work has been partially supported by the French National Research Agency (ANR) through 
COSINUS program (project COSTA-BRAVA nr. ANR-09-COSI-015). 
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